module GOLD_vert_friction
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of GOLD.                                        *
!*                                                                     *
!* GOLD is free software; you can redistribute it and/or modify it and *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* GOLD is distributed in the hope that it will be useful, but WITHOUT *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, April 1994 - October 2006                      *
!*                                                                     *
!*    This file contains the subroutine that implements vertical       *
!*  viscosity (vertvisc).                                              *
!*                                                                     *
!*    The vertical diffusion of momentum is fully implicit.  This is   *
!*  necessary to allow for vanishingly small layers.  The coupling     *
!*  is based on the distance between the centers of adjacent layers,   *
!*  except where a layer is close to the bottom compared with a        *
!*  bottom boundary layer thickness when a bottom drag law is used.    *
!*  A stress top b.c. and a no slip bottom  b.c. are used.  There      *
!*  is no limit on the time step for vertvisc.                         *
!*                                                                     *
!*    Near the bottom, the horizontal thickness interpolation scheme   *
!*  changes to an upwind biased estimate to control the effect of      *
!*  spurious Montgomery potential gradients at the bottom where        *
!*  nearly massless layers layers ride over the topography.  Within a  *
!*  few boundary layer depths of the bottom, the harmonic mean         *
!*  thickness (i.e. (2 h+ h-) / (h+ + h-) ) is used if the velocity    *
!*  is from the thinner side and the arithmetic mean thickness         *
!*  (i.e. (h+ + h-)/2) is used if the velocity is from the thicker     *
!*  side.  Both of these thickness estimates are second order          *
!*  accurate.  Above this the arithmetic mean thickness is used.       *
!*                                                                     *
!*    In addition, vertvisc truncates any velocity component that      *
!*  exceeds maxvel to truncvel. This basically keeps instabilities     *
!*  spatially localized.  The number of times the velocity is          *
!*  truncated is reported each time the energies are saved, and if     *
!*  exceeds CS%Maxtrunc the model will stop itself and change the time *
!*  to a large value.  This has proven very useful in (1) diagnosing   *
!*  model failures and (2) letting the model settle down to a          *
!*  meaningful integration from a poorly specified initial condition.  *
!*                                                                     *
!*    The same code is used for the two velocity components, by        *
!*  indirectly referencing the velocities and defining a handful of    *
!*  direction-specific defined variables.                              *
!*                                                                     *
!*  Macros written all in capital letters are defined in GOLD_memory.h.*
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v, frhatv, tauy                          *
!*    j    x ^ x ^ x   At >:  u, frhatu, taux                          *
!*    j    > o > o >   At o:  h                                        *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use GOLD_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use GOLD_diag_mediator, only : diag_ptrs, time_type
use GOLD_checksums, only : uchksum, vchksum
use GOLD_error_handler, only : GOLD_error, FATAL, WARNING
use GOLD_file_parser, only : get_param, log_version, param_file_type
use GOLD_grid, only : ocean_grid_type
use GOLD_PointAccel, only : write_u_accel, write_v_accel, PointAccel_init
use GOLD_PointAccel, only : PointAccel_CS
use GOLD_variables, only : forcing, thermo_var_ptrs, directories, vertvisc_type
use GOLD_variables, only : ocean_internal_state, ocean_OBC_type, OBC_SIMPLE

implicit none ; private

#include <GOLD_memory.h>

public vertvisc, vertvisc_remnant, vertvisc_coef
public vertvisc_limit_vel, vertvisc_init, vertvisc_end

type, public :: vertvisc_CS ; private
  real    :: Hmix           ! The mixed layer thickness in m.
  real    :: Hmix_stress    ! The mixed layer thickness over which the wind
                            ! stress is applied with direct_stress, in m.
  real    :: Kvml           ! The mixed layer vertical viscosity in m2 s-1.
  real    :: Kv             ! The interior vertical viscosity in m2 s-1.
  real    :: Hbbl           ! The static bottom boundary layer thickness, in m.
  real    :: Kvbbl          ! The vertical viscosity in the bottom boundary
                            ! layer, in m2 s-1.

  real    :: maxvel         ! Velocity components greater than maxvel,
                            ! in m s-1, are truncated.
  logical :: CFL_based_trunc ! If true, base truncations on CFL numbers, not
                            ! absolute velocities.
  real    :: CFL_trunc      ! Velocity components will be truncated when they
                            ! are large enough that the corresponding CFL number
                            ! exceeds this value, nondim. 
  real    :: CFL_report     ! The value of the CFL number that will cause the
                            ! accelerations to be reported, nondim.  CFL_report
                            ! will often equal CFL_trunc.

  real PTR_, dimension(NXMEMQP_,NYMEM_,NZp1_) :: &
    a_u                ! The u-drag coefficient across an interface, in m s-1.
  real PTR_, dimension(NXMEMQP_,NYMEM_,NZ_) :: &
    h_u                ! The effective layer thickness at u-points, m or kg m-2.
  real PTR_, dimension(NXMEM_,NYMEMQP_,NZp1_) :: &
    a_v                ! The v-drag coefficient across an interface, in m s-1.
  real PTR_, dimension(NXMEM_,NYMEMQP_,NZ_) :: &
    h_v                ! The effective layer thickness at v-points, m or kg m-2.
  real, pointer, dimension(:,:) :: &
    a1_shelf_u => NULL(), & ! The surface coupling coefficients under ice
    a1_shelf_v => NULL()    ! shelves, in m s-1. These are retained to determine
                            ! the stress under shelves.

  logical :: split          ! If true, use the split time stepping scheme.
  logical :: bottomdraglaw  ! If true, the  bottom stress is calculated with a
                            ! drag law c_drag*|u|*u. The velocity magnitude
                            ! may be an assumed value or it may be based on the
                            ! actual velocity in the bottommost HBBL, depending
                            ! on whether linear_drag is true.
  logical :: Channel_drag   ! If true, the drag is exerted directly on each
                            ! layer according to what fraction of the bottom
                            ! they overlie.
  logical :: harmonic_visc  ! If true, the harmonic mean thicknesses are used
                            ! to calculate the viscous coupling between layers
                            ! except near the bottom.  Otherwise the arithmetic
                            ! mean thickness is used except near the bottom.
  logical :: direct_stress  ! If true, the wind stress is distributed over the
                            ! topmost Hmix_stress of fluid and KVML may be very small.
  logical :: bulkmixedlayer ! If true, a refined bulk mixed layer is used with
                            ! nkml sublayers.
  logical :: dynamic_viscous_ML  ! If true, use the results from a dynamic
                            ! calculation, perhaps based on a bulk Richardson
                            ! number criterion, to determine the mixed layer
                            ! thickness for viscosity.
  logical :: debug          ! If true, write verbose checksums for debugging purposes.
  integer :: nkml           ! The number of layers in the mixed layer.
  integer, pointer :: ntrunc  ! The number of times the velocity has been
                              ! truncated since the last call to write_energy.
  character(len = 200) :: u_trunc_file ! The complete path to files in which a
  character(len = 200) :: v_trunc_file ! column's worth of accelerations are
                                       ! written when velocity truncations occur.
  type(diag_ptrs), pointer :: diag ! A pointer to a structure of shareable
                            ! ocean diagnostic fields.
  integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1
  integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1
  integer :: id_Ray_u = -1, id_Ray_v = -1, id_taux_bot = -1, id_tauy_bot = -1

  type(PointAccel_CS), pointer :: PointAccel_CSp => NULL()
end type vertvisc_CS

contains

subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, G, CS, taux_bot, tauy_bot)
!    This subroutine does a fully implicit vertical diffusion
!  of momentum.  Stress top and bottom b.c.s are used.
  real, intent(inout), dimension(NXMEMQ_,NYMEM_,NZ_) :: u
  real, intent(inout), dimension(NXMEM_,NYMEMQ_,NZ_) :: v
  real, intent(in),    dimension(NXMEM_,NYMEM_,NZ_)  :: h
  type(forcing), intent(in)                          :: fluxes
  type(vertvisc_type), intent(inout)                 :: visc
  real, intent(in)                                   :: dt
  type(ocean_OBC_type), pointer                      :: OBC
  type(ocean_grid_type), intent(in)                  :: G
  type(vertvisc_CS), pointer                         :: CS
  real, dimension(NXMEMQ_,NYMEM_), optional, intent(out) :: taux_bot
  real, dimension(NXMEM_,NYMEMQ_), optional, intent(out) :: tauy_bot
  
! Arguments: u - Zonal velocity, in m s-1.  Intent in/out.
!  (in/out)  v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      visc - The vertical viscosity type, containing information about
!                   viscosities and bottom drag-related quantities.
!  (in)      dt - Time increment in s.
!  (in)      OBC - This open boundary condition type specifies whether, where,
!                  and what open boundary conditions are used.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 vertvisc_init.
!  (out,opt) taux_bot - The zonal bottom stress from ocean to rock, in Pa.
!  (out,opt) tauy_bot - The meridional bottom stress from ocean to rock, in Pa.
!
! Fields from fluxes used in this subroutine:
!   taux: Zonal wind stress in Pa.
!   tauy: Meridional wind stress in Pa.

  real :: b1(SZIQ_(G))          ! b1 and c1 are variables used by the
  real :: c1(SZIQ_(G),SZK_(G))  ! tridiagonal solver.  c1 is nondimensional,
                                ! while b1 has units of inverse thickness.
  real :: d1(SZIQ_(G))          ! d1=1-c1 is used by the tridiagonal solver, ND.
  real :: Ray(SZIQ_(G),SZK_(G)) ! Ray is the Rayleigh-drag velocity times the
                                ! time step, in m.
  real :: b_denom_1   ! The first term in the denominator of b1, in m or kg m-2.

  real :: Hmix             ! The mixed layer thickness over which stress
                           ! is applied with direct_stress, translated into
                           ! thickness units - either m or kg m-2.
  real :: I_Hmix           ! The inverse of Hmix, in m-1 or m2 kg-1.
  real :: Idt              ! The inverse of the time step, in s-1.
  real :: dt_Rho0          ! The time step divided by the mean
                           ! density, in s m3 kg-1.
  real :: dt_m_to_H        ! The time step times the conversion from m to the
                           ! units of thickness - either s or s m3 kg-1.
  real :: h_neglect        ! A thickness that is so small it is usually lost
                           ! in roundoff and can be neglected, in m or kg m-2.

  real :: stress           !   The surface stress times the time step, divided
                           ! by the density, in units of m2 s-1.
  real :: zDS, hfr, h_a    ! Temporary variables used with direct_stress.
  real :: surface_stress(SZIQ_(G))! The same as stress, unless the wind
                           ! stress is applied as a body force, in
                           ! units of m2 s-1.

  logical :: do_i(SZIQ_(G))

  integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz
  is = G%isc ; ie = G%iec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq ; nz = G%ke

  if (.not.associated(CS)) call GOLD_error(FATAL,"GOLD_vert_friction(visc): "// &
         "Module must be initialized before it is used.")

  Hmix = CS%Hmix_stress*G%m_to_H
  I_Hmix = 1.0 / Hmix
  dt_Rho0 = dt/G%H_to_kg_m2
  dt_m_to_H = dt*G%m_to_H
  h_neglect = G%H_subroundoff
  Idt = 1.0 / dt

  do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo

  !   Update the zonal velocity component using a modification of a standard
  ! tridagonal solver.
  do j=G%jsc,G%jec
    do I=Isq,Ieq ; do_i(I) = (G%umask(I,j) > 0) ; enddo

    if (ASSOCIATED(CS%diag%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq
      CS%diag%du_dt_visc(I,j,k) = u(I,j,k)
    enddo ; enddo ; endif

!   One option is to have the wind stress applied as a body force
! over the topmost Hmix fluid.  If DIRECT_STRESS is not defined,
! the wind stress is applied as a stress boundary condition.
    if (CS%direct_stress) then
      do I=Isq,Ieq ; if (do_i(I)) then
        surface_stress(I) = 0.0
        zDS = 0.0
        stress = dt_Rho0 * fluxes%taux(I,j)
        do k=1,nz
          h_a = 0.5 * (h(I,j,k) + h(I+1,j,k)) + h_neglect
          hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a
          u(I,j,k) = u(I,j,k) + I_Hmix * hfr * stress
          zDS = zDS + h_a ; if (zDS >= Hmix) exit
        enddo
      endif ; enddo ! end of i loop
    else ; do I=Isq,Ieq
      surface_stress(I) = dt_Rho0 * (G%umask(I,j)*fluxes%taux(I,j))
    enddo ; endif ! direct_stress

    if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq
      Ray(I,k) = visc%Ray_u(I,j,k)
    enddo ; enddo ; endif

    do I=Isq,Ieq ; if (do_i(I)) then
      b_denom_1 = CS%h_u(I,j,1) + dt_m_to_H * (Ray(I,1) + CS%a_u(I,j,1))
      b1(I) = 1.0 / (b_denom_1 + dt_m_to_H*CS%a_u(I,j,2))
      d1(I) = b_denom_1 * b1(I)
      u(I,j,1) = b1(I) * (CS%h_u(I,j,1) * u(I,j,1) + surface_stress(I))
    endif ; enddo
    do k=2,nz ; do I=Isq,Ieq ; if (do_i(I)) then
      c1(I,k) = dt_m_to_H * CS%a_u(I,j,K) * b1(I)
      b_denom_1 = CS%h_u(I,j,k) + dt_m_to_H * (Ray(I,k) + CS%a_u(I,j,K)*d1(I))
      b1(I) = 1.0 / (b_denom_1 + dt_m_to_H * CS%a_u(I,j,K+1))
      d1(I) = b_denom_1 * b1(I)
      u(I,j,k) = (CS%h_u(I,j,k) * u(I,j,k) + &
                  dt_m_to_H * CS%a_u(I,j,K) * u(I,j,k-1)) * b1(I)
    endif ; enddo ; enddo
    do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then
      u(I,j,k) = u(I,j,k) + c1(I,k+1) * u(I,j,k+1)
    endif ; enddo ; enddo ! i and k loops

    if (ASSOCIATED(CS%diag%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq
      CS%diag%du_dt_visc(I,j,k) = (u(I,j,k) - CS%diag%du_dt_visc(I,j,k))*Idt
    enddo ; enddo ; endif

    if (ASSOCIATED(visc%taux_shelf)) then ; do I=Isq,Ieq
      visc%taux_shelf(I,j) = -G%Rho0*CS%a1_shelf_u(I,j)*u(I,j,1) ! - u_shelf?
    enddo ; endif
    
    if (PRESENT(taux_bot)) then
      do I=Isq,Ieq
        taux_bot(I,j) = G%Rho0 * (u(I,j,nz)*CS%a_u(I,j,nz+1))
      enddo
      if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq
        taux_bot(I,j) = taux_bot(I,j) + G%Rho0 * (Ray(I,k)*u(I,j,k))
      enddo ; enddo ; endif
    endif
  enddo ! end u-component j loop

  ! Now work on the meridional velocity component.
  do J=Jsq,Jeq
    do i=is,ie ; do_i(i) = (G%vmask(i,J) > 0) ; enddo

    if (ASSOCIATED(CS%diag%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
      CS%diag%dv_dt_visc(i,J,k) = v(i,J,k)
    enddo ; enddo ; endif

!   One option is to have the wind stress applied as a body force
! over the topmost Hmix fluid.  If DIRECT_STRESS is not defined,
! the wind stress is applied as a stress boundary condition.
    if (CS%direct_stress) then
      do i=is,ie ; if (do_i(i)) then
        surface_stress(i) = 0.0
        zDS = 0.0
        stress = dt_Rho0 * fluxes%tauy(i,J)
        do k=1,nz
          h_a = 0.5 * (h(i,J,k) + h(i,J+1,k))
          hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a
          v(i,J,k) = v(i,J,k) + I_Hmix * hfr * stress
          zDS = zDS + h_a ; if (zDS >= Hmix) exit
        enddo
      endif ; enddo ! end of i loop
    else ; do i=is,ie
      surface_stress(i) = dt_Rho0 * (G%vmask(i,J)*fluxes%tauy(i,J))
    enddo ; endif ! direct_stress

    if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie
      Ray(i,k) = visc%Ray_v(i,J,k)
    enddo ; enddo ; endif

    do i=is,ie ; if (do_i(i)) then
      b_denom_1 = CS%h_v(i,J,1) + dt_m_to_H * (Ray(i,1) + CS%a_v(i,J,1))
      b1(i) = 1.0 / (b_denom_1 + dt_m_to_H*CS%a_v(i,J,2))
      d1(i) = b_denom_1 * b1(i)
      v(i,J,1) = b1(i) * (CS%h_v(i,J,1) * v(i,J,1) + surface_stress(i))
    endif ; enddo
    do k=2,nz ; do i=is,ie ; if (do_i(i)) then
      c1(i,k) = dt_m_to_H * CS%a_v(i,J,K) * b1(i)
      b_denom_1 = CS%h_v(i,J,k) + dt_m_to_H *  (Ray(i,k) + CS%a_v(i,J,K)*d1(i))
      b1(i) = 1.0 / (b_denom_1 + dt_m_to_H * CS%a_v(i,J,K+1))
      d1(i) = b_denom_1 * b1(i)
      v(i,J,k) = (CS%h_v(i,J,k) * v(i,J,k) + dt_m_to_H * &
                  CS%a_v(i,J,K) * v(i,J,k-1)) * b1(i)
    endif ; enddo ; enddo
    do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
      v(i,J,k) = v(i,J,k) + c1(i,k+1) * v(i,J,k+1)
    endif ; enddo ; enddo ! i and k loops

    if (ASSOCIATED(CS%diag%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
      CS%diag%dv_dt_visc(i,J,k) = (v(i,J,k) - CS%diag%dv_dt_visc(i,J,k))*Idt
    enddo ; enddo ; endif

    if (ASSOCIATED(visc%tauy_shelf)) then ; do i=is,ie
      visc%tauy_shelf(i,J) = -G%Rho0*CS%a1_shelf_v(i,J)*v(i,J,1) ! - v_shelf?
    enddo ; endif

    if (present(tauy_bot)) then
      do i=is,ie
        tauy_bot(i,J) = G%Rho0 * (v(i,J,nz)*CS%a_v(i,J,nz+1))
      enddo
      if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie
        tauy_bot(i,J) = tauy_bot(i,J) + G%Rho0 * (Ray(i,k)*v(i,J,k))
      enddo ; enddo ; endif
    endif
  enddo ! end of v-component J loop

  call vertvisc_limit_vel(u, v, h, fluxes, visc, dt, G, CS)

  ! Here the velocities associated with open boundary conditions are applied.
  if (associated(OBC)) then
    if (OBC%apply_OBC_u) then
      do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq
        if (OBC%OBC_mask_u(I,j) .and. (OBC%OBC_kind_u(I,j) == OBC_SIMPLE)) &
          u(I,j,k) = OBC%u(I,j,k)
      enddo ; enddo ; enddo
    endif
    if (OBC%apply_OBC_v) then
      do k=1,nz ; do J=Jsq,Jeq ; do i=G%isc,G%iec
        if (OBC%OBC_mask_v(i,J) .and. (OBC%OBC_kind_v(i,J) == OBC_SIMPLE)) &
          v(i,J,k) = OBC%v(i,J,k)
      enddo ; enddo ; enddo
    endif
  endif
! Offer diagnostic fields for averaging.
  if (CS%id_du_dt_visc > 0) &
    call post_data(CS%id_du_dt_visc, CS%diag%du_dt_visc, CS%diag)
  if (CS%id_dv_dt_visc > 0) &
    call post_data(CS%id_dv_dt_visc, CS%diag%dv_dt_visc, CS%diag)
  if (present(taux_bot) .and. (CS%id_taux_bot > 0)) &
    call post_data(CS%id_taux_bot, taux_bot, CS%diag) 
  if (present(tauy_bot) .and. (CS%id_tauy_bot > 0)) &
    call post_data(CS%id_tauy_bot, tauy_bot, CS%diag) 

end subroutine vertvisc

subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, CS)
!    This subroutine does a fully implicit vertical diffusion
!  of momentum.  Stress top and bottom b.c.s are used.
  type(vertvisc_type), intent(in)                    :: visc
  real, intent(inout), dimension(NXMEMQ_,NYMEM_,NZ_) :: visc_rem_u
  real, intent(inout), dimension(NXMEM_,NYMEMQ_,NZ_) :: visc_rem_v
  real, intent(in)                                   :: dt
  type(ocean_grid_type), intent(in)                  :: G
  type(vertvisc_CS), pointer                         :: CS
! Arguments: visc - The vertical viscosity type, containing information about
!                   viscosities and bottom drag-related quantities, intent in.
!  (out)     visc_rem_u - Both the fraction of the momentum originally in a
!  (out)     visc_rem_v - layer that remains after a time-step of viscosity,
!                         and the fraction of a time-step's worth of a
!                         barotropic acceleration that a layer experiences
!                         after viscosity is applied, in the zonal (_u) and
!                         meridional (_v) directions.  Nondimensional between
!                         0 (at the bottom) and 1 (far above the bottom).
!  (in)      dt - Time increment in s.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 vertvisc_init.
!

  real :: b1(SZIQ_(G))          ! b1 and c1 are variables used by the
  real :: c1(SZIQ_(G),SZK_(G))  ! tridiagonal solver.  c1 is nondimensional,
                                ! while b1 has units of inverse thickness.
  real :: d1(SZIQ_(G))          ! d1=1-c1 is used by the tridiagonal solver, ND.
  real :: Ray(SZIQ_(G),SZK_(G)) ! Ray is the Rayleigh-drag velocity times the
                                ! time step, in m.
  real :: b_denom_1   ! The first term in the denominator of b1, in m or kg m-2.
  real :: dt_m_to_H        ! The time step times the conversion from m to the
                           ! units of thickness - either s or s m3 kg-1.
  logical :: do_i(SZIQ_(G))

  integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz
  is = G%isc ; ie = G%iec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq ; nz = G%ke

  if (.not.associated(CS)) call GOLD_error(FATAL,"GOLD_vert_friction(visc): "// &
         "Module must be initialized before it is used.")

  dt_m_to_H = dt*G%m_to_H

  do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo

  ! Find the zonal viscous using a modification of a standard tridagonal solver.
  do j=G%jsc,G%jec
    do I=Isq,Ieq ; do_i(I) = (G%umask(I,j) > 0) ; enddo

    if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq
      Ray(I,k) = visc%Ray_u(I,j,k)
    enddo ; enddo ; endif

    do I=Isq,Ieq ; if (do_i(I)) then
      b_denom_1 = CS%h_u(I,j,1) + dt_m_to_H * (Ray(I,1) + CS%a_u(I,j,1))
      b1(I) = 1.0 / (b_denom_1 + dt_m_to_H*CS%a_u(I,j,2))
      d1(I) = b_denom_1 * b1(I)
      visc_rem_u(I,j,1) = b1(I) * CS%h_u(I,j,1)
    endif ; enddo
    do k=2,nz ; do I=Isq,Ieq ; if (do_i(I)) then
      c1(I,k) = dt_m_to_H * CS%a_u(I,j,K)*b1(I)
      b_denom_1 = CS%h_u(I,j,k) + dt_m_to_H * (Ray(I,k) + CS%a_u(I,j,K)*d1(I))
      b1(I) = 1.0 / (b_denom_1 + dt_m_to_H * CS%a_u(I,j,K+1))
      d1(I) = b_denom_1 * b1(I)
      visc_rem_u(I,j,k) = (CS%h_u(I,j,k) + dt_m_to_H * CS%a_u(I,j,K) * visc_rem_u(I,j,k-1)) * b1(I)
    endif ; enddo ; enddo
    do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then
      visc_rem_u(I,j,k) = visc_rem_u(I,j,k) + c1(I,k+1)*visc_rem_u(I,j,k+1)
    endif ; enddo ; enddo ! i and k loops

  enddo ! end u-component j loop

  ! Now find the meridional viscous using a modification.
  do J=Jsq,Jeq
    do i=is,ie ; do_i(i) = (G%vmask(i,J) > 0) ; enddo

    if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie
      Ray(i,k) = visc%Ray_v(i,J,k)
    enddo ; enddo ; endif

    do i=is,ie ; if (do_i(i)) then
      b_denom_1 = CS%h_v(i,J,1) + dt_m_to_H * (Ray(i,1) + CS%a_v(i,J,1))
      b1(i) = 1.0 / (b_denom_1 + dt_m_to_H*CS%a_v(i,J,2))
      d1(i) = b_denom_1 * b1(i)
      visc_rem_v(i,J,1) = b1(i) * CS%h_v(i,J,1)
    endif ; enddo
    do k=2,nz ; do i=is,ie ; if (do_i(i)) then
      c1(i,k) = dt_m_to_H * CS%a_v(i,J,K)*b1(i)
      b_denom_1 = CS%h_v(i,J,k) + dt_m_to_H *  (Ray(i,k) + CS%a_v(i,J,K)*d1(i))
      b1(i) = 1.0 / (b_denom_1 + dt_m_to_H * CS%a_v(i,J,K+1))
      d1(i) = b_denom_1 * b1(i)
      visc_rem_v(i,J,k) = (CS%h_v(i,J,k) + dt_m_to_H * CS%a_v(i,J,K) * visc_rem_v(i,J,k-1)) * b1(i)
    endif ; enddo ; enddo
    do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
      visc_rem_v(i,J,k) = visc_rem_v(i,J,k) + c1(i,k+1)*visc_rem_v(i,J,k+1)
    endif ; enddo ; enddo ! i and k loops
  enddo ! end of v-component J loop


  if (CS%debug) then
    call uchksum(visc_rem_u,"visc_rem_u",G,haloshift=0)
    call vchksum(visc_rem_v,"visc_rem_v",G,haloshift=0)
  endif

end subroutine vertvisc_remnant


subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, CS)
!    This subroutine calculates the coupling coefficients (CS%a_u and CS%a_v)
! and effective layer thicknesses (CS%h_u and CS%h_v) for later use in the
! applying the implicit vertical viscosity via vertvisc.
  real, intent(in),    dimension(NXMEMQ_,NYMEM_,NZ_) :: u
  real, intent(in),    dimension(NXMEM_,NYMEMQ_,NZ_) :: v
  real, intent(in),    dimension(NXMEM_,NYMEM_,NZ_)  :: h
  type(forcing), intent(in)                          :: fluxes
  type(vertvisc_type), intent(in)                    :: visc
  real, intent(in)                                   :: dt
  type(ocean_grid_type), intent(in)                  :: G
  type(vertvisc_CS), pointer                         :: CS

! Arguments: u - Zonal velocity, in m s-1.  Intent in.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      visc - The vertical viscosity type, containing information about
!                   viscosities and bottom drag-related quantities.
!  (in)      dt - Time increment in s.
!  (in)      G - The ocean's grid structure.
!  (in/out)  CS - The control structure returned by a previous call to
!                 vertvisc_init.
!
! Field from fluxes used in this subroutine:
!   ustar: the friction velocity in m s-1, used here as the mixing
!     velocity in the mixed layer if G%nkml > 1 in a bulk mixed layer.
!
  real, dimension(SZIQ_(G),SZK_(G)) :: &
    h_harm, &   ! Harmonic mean of the thicknesses around a velocity grid point,
                ! given by 2*(h+ * h-)/(h+ + h-), in m or kg m-2 (H for short).
    hvel, &     ! hvel is the thickness used at a velocity grid point, in H.
    hvel_shelf  ! The equivalent of hvel under shelves, in H.
  real, dimension(SZIQ_(G),SZK_(G)+1) :: &
    a, &        ! The drag coefficients across interfaces, in m s-1.  a times
                ! the velocity difference gives the stress across an interface.
    a_shelf, &  ! The drag coefficients across interfaces in water columns under
                ! ice shelves, in m s-1.
    z_i         ! An estimate of each interface's height above the bottom,
                ! normalized by the bottom boundary layer thickness, nondim.
  real, dimension(SZIQ_(G)) :: &
    kv_bbl, &     ! The bottom boundary layer viscosity in m2 s-1.
    bbl_thick, &  ! The bottom boundary layer thickness in m or kg m-2.
    I_Hbbl, &     ! The inverse of the bottom boundary layer thickness, in units
                  ! of H-1 (i.e., m-1 or m2 kg-1).
    I_Htbl, &     ! The inverse of the top boundary layer thickness, in units
                  ! of H-1 (i.e., m-1 or m2 kg-1).
    zcol1, &      ! The height of the interfaces to the north and south of a
    zcol2, &      ! v-point, in m or kg m-2.
    Ztop_min, &   ! The deeper of the two adjacent surface heights, in H. 
    Dmin, &       ! The shallower of the two adjacent bottom depths converted to
                  ! thickness units, in m or kg m-2.
    zh            ! An estimate of the interface's distance from the bottom
                  ! based on harmonic mean thicknesses, in m or kg m-2.
  real, allocatable, dimension(:,:) :: hML_u, hML_v
  real :: zcol(SZI_(G)) ! The height of an interface at h-points, in m or kg m-2.
  real :: I_Hmix  ! The inverse of the mixed layer thickness, in m-1 or m2 kg-1.
  real :: botfn   ! A function which goes from 1 at the bottom to 0 much more
                  ! than Hbbl into the interior.
  real :: topfn   ! A function which goes from 1 at the top to 0 much more
                  ! than Htbl into the interior.
  real :: z2      ! The distance from the bottom, normalized by Hbbl, nondim.
  real :: h_neglect     ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected, in H.

  real :: h_arith ! The arithmetic mean thickness, in m or kg m-2.
  logical, dimension(SZIQ_(G)) :: do_i, do_i_shelf
  logical :: do_any_shelf

  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq ; nz = G%ke

  if (.not.associated(CS)) call GOLD_error(FATAL,"GOLD_vert_friction(coef): "// &
         "Module must be initialized before it is used.")

  h_neglect = G%H_subroundoff
  I_Hmix = 1.0 / (CS%Hmix * G%m_to_H + h_neglect)
  I_Hbbl(:) = 1.0 / (CS%Hbbl * G%m_to_H + h_neglect)

  if (CS%debug .or. (CS%id_hML_u > 0)) then
    allocate(hML_u(G%Isdq:G%Iedq,G%jsd:G%jed)) ; hML_u(:,:) = 0.0
  endif
  if (CS%debug .or. (CS%id_hML_v > 0)) then
    allocate(hML_v(G%isd:G%ied,G%Jsdq:G%Jedq)) ; hML_v(:,:) = 0.0
  endif

  if ((associated(visc%taux_shelf) .or. associated(fluxes%frac_shelf_u)) .and. &
      .not.associated(CS%a1_shelf_u)) then
    allocate(CS%a1_shelf_u(G%Isdq:G%Iedq,G%jsd:G%jed)) ; CS%a1_shelf_u(:,:)=0.0
  endif
  if ((associated(visc%tauy_shelf) .or. associated(fluxes%frac_shelf_u)) .and. &
      .not.associated(CS%a1_shelf_v)) then
    allocate(CS%a1_shelf_v(G%isd:G%ied,G%Jsdq:G%Jedq)) ; CS%a1_shelf_v(:,:)=0.0
  endif

  do j=G%Jsc,G%Jec
    do I=Isq,Ieq ; do_i(I) = (G%umask(I,j) > 0) ; enddo

    if (CS%bottomdraglaw) then ; do I=Isq,Ieq
      kv_bbl(I) = visc%kv_bbl_u(I,j)
      bbl_thick(I) = visc%bbl_thick_u(I,j) * G%m_to_H
      if (do_i(I)) I_Hbbl(I) = 1.0 / (bbl_thick(I) + h_neglect)
    enddo ; endif

!    The following block calculates the thicknesses at velocity
!  grid points for the vertical viscosity (hvel[k]).  Near the
!  bottom an upwind biased thickness is used to control the effect
!  of spurious Montgomery potential gradients at the bottom where
!  nearly massless layers layers ride over the topography.
    if (CS%harmonic_visc) then
      do I=Isq,Ieq ; z_i(I,nz+1) = 0.0 ; enddo
      do k=nz,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then
        h_harm(I,k) = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect)
        h_arith = 0.5*(h(i+1,j,k)+h(i,j,k))

        hvel(I,k) = h_harm(I,k)
        if (u(I,j,k) * (h(i+1,j,k)-h(i,j,k)) < 0) then
          z2 = z_i(I,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
          hvel(I,k) = (1.0-botfn)*h_harm(I,k) + botfn*h_arith
        endif
        z_i(I,k) =  z_i(I,k+1) + h_harm(I,k)*I_Hbbl(I)
      endif ; enddo ; enddo ! i & k loops
    else ! Not harmonic_visc
      do I=Isq,Ieq
        zh(I) = 0.0 ; z_i(I,nz+1) = 0.0
        Dmin(I) = min(G%D(i,j), G%D(i+1,j)) * G%m_to_H
      enddo
      do i=Isq,Ieq+1 ; zcol(i) = -G%D(i,j) * G%m_to_H ; enddo
      do k=nz,1,-1
        do i=Isq,Ieq+1 ; zcol(i) = zcol(i) + h(i,j,k) ; enddo
        do I=Isq,Ieq ; if (do_i(I)) then
          h_harm(I,k) = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect)
          h_arith = 0.5*(h(i+1,j,k)+h(i,j,k))
          zh(I) = zh(I) + h_harm(I,k)
          z_i(I,k) = max(zh(I), max(zcol(i),zcol(i+1)) + Dmin(I)) * I_Hbbl(I)

          hvel(I,k) = h_arith
          if (u(I,j,k) * (h(i+1,j,k)-h(i,j,k)) > 0) then
            z2 = max(zh(I), max(zcol(i),zcol(i+1)) + Dmin(I)) * I_Hbbl(I)
            botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
            hvel(I,k) = (1.0-botfn)*h_arith + botfn*h_harm(I,k)
          endif

        endif ; enddo ! i  loop
      enddo ! k loop
    endif

    call find_coupling_coef(a, hvel, do_i, work_on_u=.true.)
    do_any_shelf = .false.
    if (associated(fluxes%frac_shelf_u)) then
      do I=Isq,Ieq
        CS%a1_shelf_u(I,j) = 0.0
        do_i_shelf(I) = (do_i(I) .and. fluxes%frac_shelf_u(I,j) > 0.0)
        if (do_i_shelf(I)) do_any_shelf = .true.
      enddo
      if (do_any_shelf) then
        if (CS%harmonic_visc) then
          call find_coupling_coef(a_shelf, hvel, do_i_shelf, work_on_u=.true., shelf=.true.)
        else  ! Find upwind-biased thickness near the surface.
          ! Perhaps this needs to be done more carefully, via find_eta.
          do I=Isq,Ieq ; if (do_i_shelf(I)) then
            zh(I) = 0.0 ; Ztop_min(I) = min(zcol(i), zcol(i+1))
            I_HTbl(I) = 1.0 / (visc%tbl_thick_shelf_u(I,j)*G%m_to_H + h_neglect)
          endif ; enddo
          do k=1,nz
            do i=Isq,Ieq+1 ; zcol(i) = zcol(i) - h(i,j,k) ; enddo
            do I=Isq,Ieq ; if (do_i_shelf(I)) then
              h_arith = 0.5*(h(i+1,j,k)+h(i,j,k))
              zh(I) = zh(I) + h_harm(I,k)

              hvel_shelf(I,k) = hvel(I,k)
              if (u(I,j,k) * (h(i+1,j,k)-h(i,j,k)) > 0) then
                z2 = max(zh(I), Ztop_min(I) - min(zcol(i),zcol(i+1))) * I_HTbl(I)
                topfn = 1.0 / (1.0 + 0.09*z2**6)
                hvel_shelf(I,k) = min(hvel(I,k), (1.0-topfn)*h_arith + topfn*h_harm(I,k))
              endif
            endif ; enddo
          enddo
          call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, work_on_u=.true., shelf=.true.)
        endif
        do I=Isq,Ieq ; if (do_i_shelf(I)) CS%a1_shelf_u(I,j) = a_shelf(I,1) ; enddo
      endif
    endif

    if (do_any_shelf) then
      do K=1,nz+1 ; do I=Isq,Ieq ; if (do_i_shelf(I)) then
        CS%a_u(I,j,K) = fluxes%frac_shelf_u(I,j)  * a_shelf(I,k) + &
                   (1.0-fluxes%frac_shelf_u(I,j)) * a(I,K)
      elseif (do_i(I)) then
        CS%a_u(I,j,K) = a(I,K)
      endif ; enddo ; enddo
      do k=1,nz ; do I=Isq,Ieq ; if (do_i_shelf(I)) then
        ! Should we instead take the inverse of the average of the inverses?
        CS%h_u(I,j,k) = fluxes%frac_shelf_u(I,j)  * hvel_shelf(I,k) + &
                   (1.0-fluxes%frac_shelf_u(I,j)) * hvel(I,k)
      elseif (do_i(I)) then
        CS%h_u(I,j,k) = hvel(I,k)
      endif ; enddo ; enddo
    else
      do K=1,nz+1 ; do I=Isq,Ieq ; if (do_i(I)) CS%a_u(I,j,K) = a(I,K) ; enddo ; enddo
      do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) CS%h_u(I,j,k) = hvel(I,k) ; enddo ; enddo
    endif

  enddo

  ! Now work on v-points.
  do J=Jsq,Jeq
    do i=is,ie ; do_i(i) = (G%vmask(i,J) > 0) ; enddo

    if (CS%bottomdraglaw) then ; do i=is,ie
      kv_bbl(i) = visc%kv_bbl_v(i,J)
      bbl_thick(i) = visc%bbl_thick_v(i,J) * G%m_to_H
      if (do_i(i)) I_Hbbl(i) = 1.0 / bbl_thick(i)
    enddo ; endif

!    The following block calculates the thicknesses at velocity
!  grid points for the vertical viscosity (hvel[k]).  Near the
!  bottom an upwind biased thickness is used to control the effect
!  of spurious Montgomery potential gradients at the bottom where
!  nearly massless layers layers ride over the topography.
    if (CS%harmonic_visc) then
      do i=is,ie ; z_i(i,nz+1) = 0.0 ; enddo

      do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
        h_harm(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / (h(i,j,k)+h(i,j+1,k)+h_neglect)
        h_arith = 0.5*(h(i,j+1,k)+h(i,j,k))

        hvel(i,k) = h_harm(i,k)
        if (v(i,J,k) * (h(i,j+1,k)-h(i,j,k)) < 0) then
          z2 = z_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
          hvel(i,k) = (1.0-botfn)*h_harm(i,k) + botfn*h_arith
        endif
        z_i(i,k) = z_i(i,k+1)  + h_harm(i,k)*I_Hbbl(i)
      endif ; enddo ; enddo ! i & k loops
    else ! Not harmonic_visc
      do i=is,ie
        zh(i) = 0.0 ; z_i(i,nz+1) = 0.0
        Dmin(i) = min(G%D(i,j), G%D(i,j+1)) * G%m_to_H
        zcol1(i) = -G%D(i,j) * G%m_to_H ; zcol2(i) = -G%D(i,j+1) * G%m_to_H
      enddo
      do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
        h_harm(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / (h(i,j,k)+h(i,j+1,k)+h_neglect)
        h_arith = 0.5*(h(i,j+1,k)+h(i,j,k))
        zh(i) = zh(i) + h_harm(i,k)
        zcol1(i) = zcol1(i) + h(i,j,k) ; zcol2(i) = zcol2(i) + h(i,j+1,k)
        z_i(I,k) = max(zh(i), max(zcol1(i),zcol2(i)) + Dmin(i)) * I_Hbbl(i)

        hvel(i,k) = h_arith
        if (v(i,J,k) * (h(i,j+1,k)-h(i,j,k)) > 0) then
          z2 = max(zh(i), max(zcol1(i),zcol2(i)) + Dmin(i)) * I_Hbbl(i)
          botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
          hvel(i,k) = (1.0-botfn)*h_arith + botfn*h_harm(i,k)
        endif

      endif ; enddo ; enddo ! i & k loops
    endif

    call find_coupling_coef(a, hvel, do_i, work_on_u=.false.)
    do_any_shelf = .false.
    if (associated(fluxes%frac_shelf_v)) then
      do i=is,ie
        CS%a1_shelf_v(i,J) = 0.0
        do_i_shelf(i) = (do_i(i) .and. fluxes%frac_shelf_v(i,J) > 0.0)
        if (do_i_shelf(I)) do_any_shelf = .true.
      enddo
      if (do_any_shelf) then
        if (CS%harmonic_visc) then
          call find_coupling_coef(a_shelf, hvel, do_i_shelf, work_on_u=.false., shelf=.true.)
        else  ! Find upwind-biased thickness near the surface.
          ! Perhaps this needs to be done more carefully, via find_eta.
          do i=is,ie ; if (do_i_shelf(i)) then
            zh(i) = 0.0 ; Ztop_min(I) = min(zcol1(i), zcol2(i))
            I_HTbl(i) = 1.0 / (visc%tbl_thick_shelf_v(i,J)*G%m_to_H + h_neglect)
          endif ; enddo
          do k=1,nz
            do i=is,ie ; if (do_i_shelf(i)) then
              zcol1(i) = zcol1(i) - h(i,j,k) ; zcol2(i) = zcol2(i) - h(i,j+1,k)
              h_arith = 0.5*(h(i,j+1,k)+h(i,j,k))
              zh(i) = zh(i) + h_harm(i,k)

              hvel_shelf(i,k) = hvel(i,k)
              if (v(i,j,k) * (h(i,j+1,k)-h(i,j,k)) > 0) then
                z2 = max(zh(i), Ztop_min(i) - min(zcol1(i),zcol2(i))) * I_HTbl(i)
                topfn = 1.0 / (1.0 + 0.09*z2**6)
                hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith + topfn*h_harm(i,k))
              endif
            endif ; enddo
          enddo
          call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, work_on_u=.false., shelf=.true.)
        endif
        do i=is,ie ; if (do_i_shelf(i)) CS%a1_shelf_v(i,J) = a_shelf(i,1) ; enddo
      endif
    endif

    if (do_any_shelf) then
      do K=1,nz+1 ; do i=is,ie ; if (do_i_shelf(i)) then
        CS%a_v(i,J,K) = fluxes%frac_shelf_v(i,J)  * a_shelf(i,k) + &
                   (1.0-fluxes%frac_shelf_v(i,J)) * a(i,K)
      elseif (do_i(i)) then
        CS%a_v(i,J,K) = a(i,K)
      endif ; enddo ; enddo
      do k=1,nz ; do i=is,ie ; if (do_i_shelf(i)) then
        ! Should we instead take the inverse of the average of the inverses?
        CS%h_v(i,J,k) = fluxes%frac_shelf_v(i,J)  * hvel_shelf(i,k) + &
                   (1.0-fluxes%frac_shelf_v(i,J)) * hvel(i,k)
      elseif (do_i(i)) then
        CS%h_v(i,J,k) = hvel(i,k)
      endif ; enddo ; enddo
    else
      do K=1,nz+1 ; do i=is,ie ; if (do_i(i)) CS%a_v(i,J,K) = a(i,K) ; enddo ; enddo
      do k=1,nz ; do i=is,ie ; if (do_i(i)) CS%h_v(i,J,k) = hvel(i,k) ; enddo ; enddo
    endif

  enddo ! end of v-point j loop

  if (CS%debug) then
    call uchksum(CS%h_u,"vertvisc_coef h_u",G,haloshift=0)
    call vchksum(CS%h_v,"vertvisc_coef h_v",G,haloshift=0)
    call uchksum(CS%a_u,"vertvisc_coef a_u",G,haloshift=0)
    call vchksum(CS%a_v,"vertvisc_coef a_v",G,haloshift=0)
    if (allocated(hML_u)) call uchksum(hML_u,"vertvisc_coef hML_u",G,haloshift=0)
    if (allocated(hML_v)) call vchksum(hML_v,"vertvisc_coef hML_v",G,haloshift=0)
  endif

! Offer diagnostic fields for averaging.
  if (CS%id_au_vv > 0) call post_data(CS%id_au_vv, CS%a_u, CS%diag)
  if (CS%id_av_vv > 0) call post_data(CS%id_av_vv, CS%a_v, CS%diag)
  if (CS%id_h_u > 0) call post_data(CS%id_h_u, CS%h_u, CS%diag)
  if (CS%id_h_v > 0) call post_data(CS%id_h_v, CS%h_v, CS%diag)
  if (CS%id_hML_u > 0) call post_data(CS%id_hML_u, hML_u, CS%diag)
  if (CS%id_hML_v > 0) call post_data(CS%id_hML_v, hML_v, CS%diag)

  if (allocated(hML_u)) deallocate(hML_u)
  if (allocated(hML_v)) deallocate(hML_v)

  contains

  subroutine find_coupling_coef(a, hvel, do_i, work_on_u, shelf)
!    This subroutine calculates the 'coupling coefficient' (a[k]) at the
!  interfaces. If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the
!  adjacent layer thicknesses are used to calculate a[k] near the bottom.
    real,    dimension(NXMEMQ_,NZp1_), intent(out) :: a
    real,    dimension(NXMEMQ_,NZ_),   intent(in)  :: hvel
    logical, dimension(NXMEMQ_),       intent(in)  :: do_i
    logical,                           intent(in)  :: work_on_u
    logical, optional,                 intent(in)  :: shelf
! Arguments: a - The coupling coefficent across interfaces, in m.  Intent out.
!  (in)      hvel - The thickness at velocity points, in H.
!  (in)      do_i - If true, determine the a for a column.
!  (in)      work_on_u - If true, u-points are being worked on, otherwise this
!                        call is for v-points.
!  (in)      shelf - If present and true, use a surface boundary condition
!                    appropriate for an ice shelf.
    real, dimension(SZIQ_(G)) :: &
      u_star, &   ! ustar at a velocity point, in m s-1.
      absf, &     ! The average of the neighboring absolute values of f, in s-1.
      h_ml, &     ! The mixed layer depth, in m or kg m-2.
      nk_visc, &  ! The (real) interface index of the base of mixed layer.
      z_t, &      ! The distance from the top, sometimes normalized
                  ! by Hmix, in m or nondimensional.
      kv_tbl, &
      tbl_thick
    real :: h_shear ! The distance over which shears occur, m or kg m-2.
    real :: r       ! A thickness to compare with Hbbl, in m or kg m-2.
    real :: visc_ml ! The mixed layer viscosity, in m2 s-1.
    real :: a_ml    ! The layer coupling coefficient across an interface in
                    ! the mixed layer, in m s-1.
    real :: temp1   ! A temporary variable in m2 s-1.
    real :: h_neglect   ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected, in H.
    real :: dz_neglect  ! A thickness in m that is so small it is usually lost
                        ! in roundoff and can be neglected, in m.
    real :: z2      ! A copy of z_i, nondim.
    real :: topfn
    real :: a_top
    logical :: do_shelf
    integer :: i, k, is, ie, max_nk

    if (work_on_u) then ; is = G%Iscq ; ie = G%Iecq
    else ; is = G%isc ; ie = G%iec ; endif
    h_neglect = G%H_subroundoff
    dz_neglect = G%H_subroundoff*G%H_to_m

    do_shelf = .false. ; if (present(shelf)) do_shelf = shelf

!    The following loop calculates the vertical average velocity and
!  surface mixed layer contributions to the vertical viscosity.
    do i=is,ie ; a(i,1) = 0.0 ; enddo
    if ((G%nkml>0) .or. do_shelf) then ; do k=2,nz ; do i=is,ie
      if (do_i(i)) a(i,K) = 2.0*CS%Kv
    enddo ; enddo ; else
      do i=is,ie ; z_t(i) = h_neglect*I_Hmix ; enddo
      do K=2,nz ; do i=is,ie ; if (do_i(i)) then
        z_t(i) = z_t(i) + h_harm(i,k-1)*I_Hmix
        a(i,K) = 2.0*CS%Kv + 2.0*CS%Kvml / ((z_t(i)*z_t(i)) *  &
                 (1.0 + 0.09*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)))
      endif ; enddo ; enddo
    endif

    do i=is,ie ; if (do_i(i)) then
      if (CS%bottomdraglaw) then
        r = hvel(i,nz)*0.5
        if (r < bbl_thick(i)) then
          a(i,nz+1) = 1.0*kv_bbl(i) / (1e-10*dt*kv_bbl(i) + r*G%H_to_m)
        else
          a(i,nz+1) = 1.0*kv_bbl(i) / (1e-10*dt*kv_bbl(i) + bbl_thick(i)*G%H_to_m)
        endif
      else
        a(i,nz+1) = 2.0*CS%Kvbbl / (hvel(i,nz)*G%H_to_m + 2.0e-10*dt*CS%Kvbbl)
      endif
    endif ; enddo

    do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then
      !    botfn determines when a point is within the influence of the bottom
      !  boundary layer, going from 1 at the bottom to 0 in the interior.
      z2 = z_i(i,k)
      botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)

      if (CS%bottomdraglaw) then
        a(i,K) = a(i,K) + 2.0*(kv_bbl(i)-CS%Kv)*botfn
        r = (hvel(i,k)+hvel(i,k-1))
        if (r > 2.0*bbl_thick(i)) then
          h_shear = ((1.0 - botfn) * r + botfn*2.0*bbl_thick(i))
        else
          h_shear = r
        endif
      else
        a(i,K) = a(i,K) + 2.0*(CS%Kvbbl-CS%Kv)*botfn
        h_shear = hvel(i,k) + hvel(i,k-1) + h_neglect
      endif
      !   The term including 1e-10 in the denominators is here to avoid
      ! truncation error problems in the tridiagonal solver. Effectively, this
      ! sets the maximum coupling coefficient at 1e10 m.
      a(i,K) = a(i,K) / (h_shear*G%H_to_m + 1.0e-10*dt*a(i,K))
    endif ; enddo ; enddo ! i & k loops

    if (do_shelf) then
      ! Set the coefficients to include the no-slip surface stress.
      do i=is,ie ; if (do_i(i)) then
        if (work_on_u) then
          kv_tbl(i) = visc%kv_tbl_shelf_u(I,j)
          tbl_thick(i) = visc%tbl_thick_shelf_u(I,j) * G%m_to_H
        else
          kv_tbl(i) = visc%kv_tbl_shelf_v(i,J)
          tbl_thick(i) = visc%tbl_thick_shelf_v(i,J) * G%m_to_H
        endif
        z_t(i) = 0.0

        ! If a(i,1) were not already 0, it would be added here.
        if (0.5*hvel(i,1) > tbl_thick(i)) then
          a(i,1) = kv_tbl(i) / (tbl_thick(i) *G%H_to_m + (1.0e-10*dt)*kv_tbl(i))
        else
          a(i,1) = kv_tbl(i) / (0.5*hvel(i,1)*G%H_to_m + (1.0e-10*dt)*kv_tbl(i))
        endif
      endif ; enddo

      do K=2,nz ; do i=is,ie ;  if (do_i(i)) then
        z_t(i) = z_t(i) + hvel(i,k-1) / tbl_thick(i)
        topfn = 1.0 / (1.0 + 0.09 * z_t(i)**6)

        r = (hvel(i,k)+hvel(i,k-1))
        if (r > 2.0*tbl_thick(i)) then
          h_shear = ((1.0 - topfn) * r + topfn*2.0*tbl_thick(i))
        else
          h_shear = r
        endif
      !   The term including 1e-10 in the denominators is here to avoid
      ! truncation error problems in the tridiagonal solver. Effectively, this
      ! sets the maximum coupling coefficient increment to 1e10 m.
        a_top = 2.0 * topfn * kv_tbl(i)
        a(i,K) = a(i,K) + a_top / (h_shear*G%H_to_m + 1.0e-10*dt*a_top)
      endif ; enddo ; enddo
    elseif (CS%dynamic_viscous_ML .or. (G%nkml>0)) then
      max_nk = 0
      do i=is,ie ; if (do_i(i)) then
        if (G%nkml>0) nk_visc(i) = real(G%nkml+1)
        if (work_on_u) then
          u_star(I) = 0.5*(fluxes%ustar(i,j) + fluxes%ustar(i+1,j))
          absf(I) = 0.5*(abs(G%f(I,J-1)) + abs(G%f(I,J)))
          if (CS%dynamic_viscous_ML) nk_visc(I) = visc%nkml_visc_u(I,j) + 1
        else
          u_star(i) = 0.5*(fluxes%ustar(i,j) + fluxes%ustar(i,j+1))
          absf(i) = 0.5*(abs(G%f(I-1,J)) + abs(G%f(I,J)))
          if (CS%dynamic_viscous_ML) nk_visc(i) = visc%nkml_visc_v(i,J) + 1
        endif
        h_ml(i) = h_neglect ; z_t(i) = 0.0
        max_nk = max(max_nk,ceiling(nk_visc(i) - 1.0))
      endif ; enddo

      do k=1,max_nk ; do i=is,ie ; if (do_i(i)) then
        if (k+1 <= nk_visc(i)) then ! This layer is all in the ML.
          h_ml(i) = h_ml(i) + hvel(i,k)
        elseif (k < nk_visc(i)) then ! Part of this layer is in the ML.
          h_ml(i) = h_ml(i) + (nk_visc(i) - k) * hvel(i,k)
        endif
      endif ; enddo ; enddo
      if (work_on_u .and. allocated(hML_u)) then
        do i=is,ie ; if (do_i(i)) then ; hML_u(I,j) = h_ml(I) ; endif ; enddo
      elseif (.not.work_on_u .and. allocated(hML_v)) then
        do i=is,ie ; if (do_i(i)) then ; hML_v(i,J) = h_ml(i) ; endif ; enddo
      endif

      do K=2,max_nk ; do i=is,ie ; if (do_i(i)) then ; if (k < nk_visc(i)) then
        ! Set the viscosity at the interfaces.
        z_t(i) = z_t(i) + hvel(i,k-1)
        temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i)) * G%H_to_m
        !   This viscosity is set to go to 0 at the mixed layer top and bottom
        ! (in a log-layer) and be further limited by rotation to give the
        ! natural Ekman length.
        visc_ml = u_star(i) * 0.41 * (temp1*u_star(i)) / &
                       (absf(i)*temp1 + h_ml(i)*u_star(i))
        a_ml = 4.0*visc_ml / ((hvel(i,k)+hvel(i,k-1) + h_neglect) * G%H_to_m + &
                              2.0e-10*dt*visc_ml)
        ! Choose the largest estimate of a.
        if (a_ml > a(i,K)) a(i,K) = a_ml
      endif ; endif ; enddo ; enddo
    endif

  end subroutine find_coupling_coef

end subroutine vertvisc_coef

subroutine vertvisc_limit_vel(u, v, h, fluxes, visc, dt, G, CS)
!  Within this subroutine, velocity components which exceed a threshold for
! physically reasonable values are truncated. Optionally, any column with
! excessive velocities may be sent to a diagnostic reporting subroutine.
  real, intent(inout), dimension(NXMEMQ_,NYMEM_,NZ_) :: u
  real, intent(inout), dimension(NXMEM_,NYMEMQ_,NZ_) :: v
  real, intent(in),    dimension(NXMEM_,NYMEM_,NZ_)  :: h
  type(forcing), intent(in)                          :: fluxes
  type(vertvisc_type), intent(in)                    :: visc
  real, intent(in)                                   :: dt
  type(ocean_grid_type), intent(in)                  :: G
  type(vertvisc_CS), pointer                         :: CS

  real :: maxvel           ! Velocities components greater than maxvel
  real :: truncvel         ! are truncated to truncvel, both in m s-1.
  real :: CFL              ! The local CFL number.
  real :: vel_report(SZIQ_(G))
  logical :: trunc_any, dowrite(SZIQ_(G))
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq

  maxvel = CS%maxvel
  truncvel = 0.9*maxvel

  if (len_trim(CS%u_trunc_file) > 0) then
    if (.not.CS%CFL_based_trunc) vel_report(:) = maxvel

    do j=js,je
      trunc_any = .false.
      do I=Isq,Ieq ; dowrite(I) = .false. ; enddo
      if (CS%CFL_based_trunc) then
        do I=Isq,Ieq ; vel_report(:) = 3.0e8 ; enddo ! Speed of light default.
        do k=1,nz ; do I=Isq,Ieq
          if (u(I,j,k) < 0.0) then
            CFL = (-u(I,j,k) * dt) * (G%dy_u(I,j) * G%IDXDYh(i+1,j))
          else
            CFL = (u(I,j,k) * dt) * (G%dy_u(I,j) * G%IDXDYh(i,j))
          endif
          if (CFL > CS%CFL_trunc) trunc_any = .true.
          if (CFL > CS%CFL_report) then
            dowrite(I) = .true.
            vel_report(I) = MIN(vel_report(I), abs(u(I,j,k)))
          endif
        enddo ; enddo
      else
        do k=1,nz ; do I=Isq,Ieq ; if (abs(u(I,j,k)) > maxvel) then
          dowrite(I) = .true. ; trunc_any = .true.
        endif ;enddo ; enddo
      endif
      
      do I=Isq,Ieq ; if (dowrite(I)) then
!   Here the diagnostic reporting subroutines are called if
! unphysically large values were found.
        call write_u_accel(I, j, u, h, dt, G, CS%PointAccel_CSp, &
               vel_report(I), -vel_report(I), (dt*fluxes%taux(I,j)/G%Rho0), &
               a=CS%a_u(:,j,:), hv=CS%h_u(:,j,:))
      endif ; enddo
      if (trunc_any) then ; if (CS%CFL_based_trunc) then
        do k=1,nz ; do I=Isq,Ieq
          if ((u(I,j,k) * (dt * G%dy_u(I,j))) * G%IDXDYh(i+1,j) < -CS%CFL_trunc) then
            u(I,j,k) = (-0.9*CS%CFL_trunc) * (G%DXDYh(i+1,j) / (dt * G%dy_u(I,j)))
            if (h(i,j,k) + h(i+1,j,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
          elseif ((u(I,j,k) * (dt * G%dy_u(I,j))) * G%IDXDYh(i,j) > CS%CFL_trunc) then
            u(I,j,k) = (0.9*CS%CFL_trunc) * (G%DXDYh(i,j) / (dt * G%dy_u(I,j)))
            if (h(i,j,k) + h(i+1,j,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
          endif
        enddo ; enddo
      else
        do k=1,nz ; do I=Isq,Ieq ; if (abs(u(I,j,k)) > maxvel) then
          u(I,j,k) = SIGN(truncvel,u(I,j,k))
          if (h(i,j,k) + h(i+1,j,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
        endif ; enddo ;  enddo
      endif ; endif
    enddo ! j-loop
  else
    if (CS%CFL_based_trunc) then
      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
        if ((u(I,j,k) * (dt * G%dy_u(I,j))) * G%IDXDYh(i+1,j) < -CS%CFL_trunc) then
          u(I,j,k) = (-0.9*CS%CFL_trunc) * (G%DXDYh(i+1,j) / (dt * G%dy_u(I,j)))
          if (h(i,j,k) + h(i+1,j,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
        elseif ((u(I,j,k) * (dt * G%dy_u(I,j))) * G%IDXDYh(i,j) > CS%CFL_trunc) then
          u(I,j,k) = (0.9*CS%CFL_trunc) * (G%DXDYh(i,j) / (dt * G%dy_u(I,j)))
          if (h(i,j,k) + h(i+1,j,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
        endif
      enddo ; enddo ; enddo
    else
      do k=1,nz ; do j=js,je ; do I=Isq,Ieq ; if (abs(u(I,j,k)) > maxvel) then
        u(I,j,k) = SIGN(truncvel,u(I,j,k))
        if (h(i,j,k) + h(i+1,j,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
      endif ; enddo ; enddo ; enddo
    endif
  endif

  if (len_trim(CS%v_trunc_file) > 0) then
    if (.not.CS%CFL_based_trunc) vel_report(:) = maxvel

    do J=Jsq,Jeq
      trunc_any = .false.
      do i=is,ie ; dowrite(i) = .false. ; enddo
      if (CS%CFL_based_trunc) then
        do i=is,ie ; vel_report(:) = 3.0e8 ; enddo ! Speed of light default.
        do k=1,nz ; do i=is,ie
          if (v(i,J,k) < 0.0) then
            CFL = (-v(i,J,k) * dt) * (G%dx_v(i,J) * G%IDXDYh(i,j+1))
          else
            CFL = (v(i,J,k) * dt) * (G%dx_v(i,J) * G%IDXDYh(i,j))
          endif
          if (CFL > CS%CFL_trunc) trunc_any = .true.
          if (CFL > CS%CFL_report) then
            dowrite(i) = .true.
            vel_report(i) = MIN(vel_report(i), abs(v(i,J,k)))
          endif
        enddo ; enddo
      else
        do k=1,nz ; do i=is,ie ; if (abs(v(i,J,k)) > maxvel) then
          dowrite(i) = .true. ; trunc_any = .true.
        endif ; enddo ; enddo
      endif

      do i=is,ie ; if (dowrite(i)) then
!   Here the diagnostic reporting subroutines are called if
! unphysically large values were found.
        call write_v_accel(i, J, v, h, dt, G, CS%PointAccel_CSp, &
               vel_report(I), -vel_report(I), (dt*fluxes%tauy(i,J)/G%Rho0), &
               a=CS%a_v(:,J,:),hv=CS%h_v(:,J,:))
      endif ; enddo
      if (trunc_any) then ; if (CS%CFL_based_trunc) then
        do k=1,nz; do i=is,ie
          if ((v(i,J,k) * (dt * G%dx_v(i,J))) * G%IDXDYh(i,j+1) < -CS%CFL_trunc) then
            v(i,J,k) = (-0.9*CS%CFL_trunc) * (G%DXDYh(i,j+1) / (dt * G%dx_v(i,J)))
            if (h(i,j,k) + h(i,j+1,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
          elseif ((v(i,J,k) * (dt * G%dx_v(i,J))) * G%IDXDYh(i,j) > CS%CFL_trunc) then
            v(i,J,k) = (0.9*CS%CFL_trunc) * (G%DXDYh(i,j) / (dt * G%dx_v(i,J)))
            if (h(i,j,k) + h(i,j+1,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
          endif
        enddo ; enddo
      else
        do k=1,nz ; do i=is,ie ; if (abs(v(i,J,k)) > maxvel) then
          v(i,J,k) = SIGN(truncvel,v(i,J,k))
          if (h(i,j,k) + h(i,j+1,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
        endif ; enddo ;  enddo
      endif ; endif
    enddo ! J-loop
  else
    if (CS%CFL_based_trunc) then
      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
        if ((v(i,J,k) * (dt * G%dx_v(i,J))) * G%IDXDYh(i,j+1) < -CS%CFL_trunc) then
          v(i,J,k) = (-0.9*CS%CFL_trunc) * (G%DXDYh(i,j+1) / (dt * G%dx_v(i,J)))
          if (h(i,j,k) + h(i,j+1,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
        elseif ((v(i,J,k) * (dt * G%dx_v(i,J))) * G%IDXDYh(i,j) > CS%CFL_trunc) then
          v(i,J,k) = (0.9*CS%CFL_trunc) * (G%DXDYh(i,j) / (dt * G%dx_v(i,J)))
          if (h(i,j,k) + h(i,j+1,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
        endif
      enddo ; enddo ; enddo
    else
      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (abs(v(i,J,k)) > maxvel) then
        v(i,J,k) = SIGN(truncvel,v(i,J,k))
        if (h(i,j,k) + h(i,j+1,k) > 6.0*G%Angstrom) CS%ntrunc = CS%ntrunc + 1
      endif ; enddo ; enddo ; enddo
    endif
  endif

end subroutine vertvisc_limit_vel

subroutine vertvisc_init(HIS, Time, G, param_file, diag, dirs, ntrunc, CS)
  type(ocean_internal_state), target, intent(in) :: HIS
  type(time_type), target, intent(in)    :: Time
  type(ocean_grid_type),   intent(in)    :: G
  type(param_file_type),   intent(in)    :: param_file
  type(diag_ptrs), target, intent(inout) :: diag
  type(directories),       intent(in)    :: dirs
  integer, target,         intent(inout) :: ntrunc
  type(vertvisc_CS),       pointer       :: CS
! Arguments: HIS - For "GOLD Internal State" a set of pointers to the fields and
!                  accelerations that make up the ocean's physical state.
!  (in)      Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure containing pointers to common diagnostic fields.
!  (in)      dirs - A structure containing several relevant directory paths.
!  (in/out)  ntrunc - The integer that stores the number of times the velocity
!                     has been truncated since the last call to write_energy.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module

  real :: hmix_str_dflt
  integer :: isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq, nz
  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  character(len=40)  :: mod = "GOLD_vert_friction" ! This module's name.
  character(len=40)  :: thickness_units = "meters or kg m-2"

  if (associated(CS)) then
    call GOLD_error(WARNING, "vertvisc_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(CS)

  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

  CS%diag => diag ; CS%ntrunc => ntrunc ; ntrunc = 0

! Default, read and log parameters
  call log_version(param_file, mod, version, tagname, "")
  call get_param(param_file, mod, "BOTTOMDRAGLAW", CS%bottomdraglaw, &
                 "If true, the bottom stress is calculated with a drag \n"//&
                 "law of the form c_drag*|u|*u. The velocity magnitude \n"//&
                 "may be an assumed value or it may be based on the \n"//&
                 "actual velocity in the bottommost HBBL, depending on \n"//&
                 "LINEAR_DRAG.", default=.true.)
  call get_param(param_file, mod, "CHANNEL_DRAG", CS%Channel_drag, &
                 "If true, the bottom drag is exerted directly on each \n"//&
                 "layer proportional to the fraction of the bottom it \n"//&
                 "overlies.", default=.false.)
  call get_param(param_file, mod, "DIRECT_STRESS", CS%direct_stress, &
                 "If true, the wind stress is distributed over the \n"//&
                 "topmost HMIX of fluid (like in HYCOM), and KVML may be \n"//&
                 "set to a very small value.", default=.false.)
  call get_param(param_file, mod, "DYNAMIC_VISCOUS_ML", CS%dynamic_viscous_ML, &
                 "If true, use a bulk Richardson number criterion to \n"//&
                 "determine the mixed layer thickness for viscosity.", &
                 default=.false.)
  call get_param(param_file, mod, "U_TRUNC_FILE", CS%u_trunc_file, &
                 "The absolute path to a file into which the accelerations \n"//&
                 "leading to zonal velocity truncations are written.  \n"//&
                 "Undefine this for efficiency if this diagnostic is not \n"//&
                 "needed.", default=" ")
  call get_param(param_file, mod, "V_TRUNC_FILE", CS%v_trunc_file, &
                 "The absolute path to a file into which the accelerations \n"//&
                 "leading to meridional velocity truncations are written. \n"//&
                 "Undefine this for efficiency if this diagnostic is not \n"//&
                 "needed.", default=" ")
  call get_param(param_file, mod, "HARMONIC_VISC", CS%harmonic_visc, &
                 "If true, use the harmonic mean thicknesses for \n"//&
                 "calculating the vertical viscosity.", default=.false.)
  call get_param(param_file, mod, "DEBUG", CS%debug, default=.false.)

  CS%Hmix_stress = -1.0e6
  if (G%nkml < 1) &
    call get_param(param_file, mod, "HMIX", CS%Hmix, &
                 "The depth over which the near-surface viscosity is \n"//&
                 "elevated if BULKMIXEDLAYER is false.", units="m", &
                 default=-1.0e6, fail_if_missing=.true.)
  if (CS%direct_stress) then
    call get_param(param_file, mod, "HMIX", hmix_str_dflt, default=-1.0e6)
    call get_param(param_file, mod, "HMIX_STRESS", CS%Hmix_stress, &
                 "The depth over which the wind stress is applied if \n"//&
                 "DIRECT_STRESS is true.", units="m", default=hmix_str_dflt)
    if (CS%Hmix_stress <= 0.0) call GOLD_error(FATAL, "vertvisc_init: " // &
       "HMIX_STRESS must be set to a positive value if DIRECT_STRESS is true.")
  endif
  call get_param(param_file, mod, "KV", CS%Kv, &
                 "The background kinematic viscosity in the interior. \n"//&
                 "The molecular value, ~1e-6 m2 s-1, may be used.", &
                 units="m2 s-1", fail_if_missing=.true.)

! CS%Kvml = CS%Kv ; CS%Kvbbl = CS%Kv ! Needed? -AJA
  if (G%nkml < 1) call get_param(param_file, mod, "KVML", CS%Kvml, &
                 "The kinematic viscosity in the mixed layer.  A typical \n"//&
                 "value is ~1e-2 m2 s-1. KVML is not used if \n"//&
                 "BULKMIXEDLAYER is true.  The default is set by KV.", &
                 units="m2 s-1", default=CS%Kv)
  if (.not.CS%bottomdraglaw) call get_param(param_file, mod, "KVBBL", CS%Kvbbl, &
                 "The kinematic viscosity in the benthic boundary layer. \n"//&
                 "A typical value is ~1e-2 m2 s-1. KVBBL is not used if \n"//&
                 "BOTTOMDRAGLAW is true.  The default is set by KV.", &
                 units="m2 s-1", default=CS%Kv)
  call get_param(param_file, mod, "HBBL", CS%Hbbl, &
                 "The thickness of a bottom boundary layer with a \n"//&
                 "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or \n"//&
                 "the thickness over which near-bottom velocities are \n"//&
                 "averaged for the drag law if BOTTOMDRAGLAW is defined \n"//&
                 "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true.)
  call get_param(param_file, mod, "MAXVEL", CS%maxvel, &
                 "The maximum velocity allowed before the velocity \n"//&
                 "components are truncated.", units="m s-1", default=3.0e8)
  call get_param(param_file, mod, "CFL_BASED_TRUNCATIONS", CS%CFL_based_trunc, &
                 "If true, base truncations on the CFL number, and not an \n"//&
                 "absolute speed.", default=.true.)
  call get_param(param_file, mod, "CFL_TRUNCATE", CS%CFL_trunc, &
                 "The value of the CFL number that will cause velocity \n"//&
                 "components to be truncated; instability can occur past 0.5.", &
                 units="nondim", default=0.5)
  call get_param(param_file, mod, "CFL_REPORT", CS%CFL_report, &
                 "The value of the CFL number that causes accelerations \n"//&
                 "to be reported; the default is CFL_TRUNCATE.", &
                 units="nondim", default=CS%CFL_trunc)

  ALLOC(CS%a_u(Isdq:Iedq,jsd:jed,nz+1)) ; CS%a_u(:,:,:) = 0.0
  ALLOC(CS%h_u(Isdq:Iedq,jsd:jed,nz))   ; CS%h_u(:,:,:) = 0.0
  ALLOC(CS%a_v(isd:ied,Jsdq:Jedq,nz+1)) ; CS%a_v(:,:,:) = 0.0
  ALLOC(CS%h_v(isd:ied,Jsdq:Jedq,nz))   ; CS%h_v(:,:,:) = 0.0

  CS%id_au_vv = register_diag_field('ocean_model', 'au_visc', G%axesui, Time, &
     'Zonal Viscous Vertical Coupling Coefficient', 'meter second-1')
  CS%id_av_vv = register_diag_field('ocean_model', 'av_visc', G%axesvi, Time, &
     'Meridional Viscous Vertical Coupling Coefficient', 'meter second-1')

  CS%id_h_u = register_diag_field('ocean_model', 'Hu_visc', G%axesuL, Time, &
     'Thickness at Zonal Velocity Points for Viscosity', thickness_units)
  CS%id_h_v = register_diag_field('ocean_model', 'Hv_visc', G%axesvL, Time, &
     'Thickness at Meridional Velocity Points for Viscosity', thickness_units)
  CS%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', G%axesu1, Time, &
     'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units)
  CS%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', G%axesv1, Time, &
     'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units)

  CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', G%axesul, &
     Time, 'Zonal Acceleration from Vertical Viscosity', 'meter second-2')
  if (CS%id_du_dt_visc > 0) call safe_alloc_ptr(diag%du_dt_visc,Isdq,Iedq,jsd,jed,nz)
  CS%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', G%axesvl, &
     Time, 'Meridional Acceleration from Vertical Viscosity', 'meter second-2')
  if (CS%id_dv_dt_visc > 0) call safe_alloc_ptr(diag%dv_dt_visc,isd,ied,Jsdq,Jedq,nz)

  CS%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', G%axesu1, &
     Time, 'Zonal Bottom Stress from Ocean to Earth', 'Pa')
  CS%id_tauy_bot = register_diag_field('ocean_model', 'tauy_bot', G%axesv1, &
     Time, 'Meridional Bottom Stress from Ocean to Earth', 'Pa')

  if ((len_trim(CS%v_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &
    call PointAccel_init(HIS, Time, G, param_file, diag, dirs, CS%PointAccel_CSp)

end subroutine vertvisc_init

subroutine vertvisc_end(CS)
  type(vertvisc_CS),   pointer       :: CS
  DEALLOC(CS%a_u) ; DEALLOC(CS%h_u)
  DEALLOC(CS%a_v) ; DEALLOC(CS%h_v)
  if (associated(CS%a1_shelf_u)) deallocate(CS%a1_shelf_u)
  if (associated(CS%a1_shelf_v)) deallocate(CS%a1_shelf_v)
  deallocate(CS)
end subroutine vertvisc_end

end module GOLD_vert_friction
